/*
 *  viscoelastic.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 1/04/11.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "viscoelastic.h"

Viscoelastic::Viscoelastic()
:	Material(),
	vTau(),
	vBeta(),
	nTau(),
	matHyper()
{}

// ***

Viscoelastic::Viscoelastic(ifstream &in, const string& na, const string & t, Material* ma)
:	Material(na,t),
	vTau(),
	vBeta(),
	vMu(),
	vEta(),
	nTau(),
	mode(0),
	rheological(0),
	matHyper(ma)
{
	int i,j;
	string keyword;
	string line;
	vector<string>	tauValues;
	stringstream stnum;
	double value;
	
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
                if(!keyword.compare("RHEOLOGICAL"))
                {
                        rheological = 1;
                }
                //in >> keyword;
                if(!keyword.compare("PREINTEGRATED"))
                {
                        mode = 1;
                }
                //in >> keyword
                if(!keyword.compare("TAU"))
		{
			getline(in,line);
			istringstream iss(line);
			while (iss >> keyword )
			{
				tauValues.push_back(keyword);
			}
			nTau = tauValues.size();
			for(i=0; i<nTau; ++i)	// Starts from 1 because the first element is the keyword V0
			{
				stnum << tauValues[i];
				stnum >> value;
				vTau.push_back(value);
				stnum.clear();
			}
		}
		//in >> keyword;
		if(!keyword.compare("BETA"))
		{
			for (i=0; i<nTau; ++i)
			{
				in >> value;
				vBeta.push_back(value);
			}
		}
        }
        /*      if(!keyword.compare("MU"))
                {
                        mode = 0;
                        for (i=0; i<nTau; ++i)
                        {
                                in >> value;
                                vMu.push_back(value);
                        }
                } */
	
	
        hypertype = matHyper->hypertype;
	
        vEta.setZero(nTau);
        vMu.setZero(nTau);
        
        if (mode == 0)
        {
                Tv.setZero(6*nTau,6);
                for(i=0;i<nTau;++i)
                {
                        vMu[i] = vBeta[i] * matHyper->getE() / 2.0;
                        vEta[i] = 2.0 * vMu[i] * vTau[i];
                        Tv.block(6*i,0,6,6) = - 2.0 * vMu[i] * MatrixXd::Identity(6,6);
                }
        }
        
        if (rheological == 1)
        {
                I.setZero(6);
                I.head(3) << 1.0,
                             1.0,
                             1.0;
        }
}

// ***

double Viscoelastic::getW( const VectorXd & C , const MatrixXd & Gamma )
{
        if(!matHyper->hypertype.compare("coupled"))
        { 
                W = matHyper->getW(C);

                switch (rheological)
                {
                        case 0:
                                for(int i=0; i < nTau; ++i)
                                        W += vMu[i] * dotVoigt(Gamma.col(i), Gamma.col(i)) - vBeta[i] * dotVoigt(matHyper->getS(C), Gamma.col(i)) + vBeta[i] * matHyper->getW(C);
                                break;
                        case 1:
                                for(int i=0; i < nTau; ++i)
                                        W += vMu[i] * dotVoigt(0.5 * (C - I) - Gamma.col(i), 0.5 * (C - I) - Gamma.col(i));
                                break;
                }
        }
        else
        {
                W = matHyper->getW(C);
                
                switch (rheological)
                {
                        case 0:
                                for(int i=0; i < nTau; ++i)
                                        W += vMu[i] * dotVoigt(Gamma.col(i), Gamma.col(i)) - vBeta[i] * dotVoigt(matHyper->getS_Iso(C), Gamma.col(i)) + vBeta[i] * matHyper->getW_Iso(C);
                                break;
                        case 1:
                                for(int i=0; i < nTau; ++i)
                                        W +=  vMu[i] * dotVoigt(0.5 * (matHyper->getJ(C) * C - I)- Gamma.col(i), 0.5 * (matHyper->getJ(C)*C - I)- Gamma.col(i));
                                break;
                }
        }
	
	return W;
}

VectorXd Viscoelastic::getS( const VectorXd & C, const MatrixXd & Gamma ) 
{
	if(matHyper->type.compare("svk")) matHyper->calcCt(C); //Means that hyperelastic material is not svk.
	
        if(!matHyper->hypertype.compare("coupled"))
        {
                S = matHyper->getS(C);
                
                switch (rheological)
                {
                        case 0:
                                for(int i=0;i<nTau;++i)
                                         S += vBeta[i] * ( matHyper->getS(C) - matHyper->Ct * makeVoigt(Gamma.col(i)) );
                                         //S += vBeta[i] * ( matHyper->getS(C) - matHyper->Ct * Gamma.col(i) );
                                break;
                        case 1:
                                for(int i=0;i<nTau;++i)
                                         S += getQv(i,C,Gamma);
                                break;
                }
        }
        else
        {
                S = matHyper->getS(C);
                
                switch (rheological)
                {
                        case 0:
                                for(int i=0;i<nTau;++i)
                                        S += vBeta[i] * ( matHyper->getS_Iso(C) - matHyper->Ctiso * makeVoigt(Gamma.col(i)));
                                        //S += vBeta[i] * ( matHyper->getS_Iso(C) - matHyper->Ctiso * Gamma.col(i) );
                                break;
                        case 1:
                                for(int i=0;i<nTau;++i)
                                        S += getQv(i,C,Gamma);
                                break;

                }
        }
	return S;
}

//***

VectorXd Viscoelastic::getQv(int i, const VectorXd & C, const MatrixXd & Gamma )
{
        if(!matHyper->hypertype.compare("coupled"))
        {
                switch (rheological)
                {
                        case 0:
                                return - ( 2.0 * vMu[i] * Gamma.col(i) - vBeta[i] * matHyper->getS(C) );
                                break;
                        case 1: 
                                return 2.0 * vMu[i] * (0.5 * (C - I) - Gamma.col(i));
                                break;
                }
        }
        else
        {
                switch (rheological)
                {
                        case 0:
                                return - ( 2.0 * vMu[i] * Gamma.col(i) - vBeta[i] * matHyper->getS_Iso(C) );
                                break;
                        case 1:
                                return 2.0 * vMu[i] * (0.5 * (matHyper->getJ(C) * C - I) - Gamma.col(i));
                                break;
                }
        }
}

// ***

VectorXd Viscoelastic::getS_Vol( const VectorXd & C, const MatrixXd & Gamma )
{
	return matHyper->getS_Vol(C);
}

//***

VectorXd Viscoelastic::getS_Iso( const VectorXd & C, const MatrixXd & Gamma)
{
	return matHyper->getS_Vol(C);
}

//************************************* Hyperelastic Response *********************************************

VectorXd Viscoelastic::getS( const VectorXd & C)
{
	return matHyper->getS(C);
}

// ***

VectorXd Viscoelastic::getS_Vol( const VectorXd & C)
{
	return matHyper->getS_Vol(C);
}

//***

VectorXd Viscoelastic::getS_Iso( const VectorXd & C)
{
	return matHyper->getS_Iso(C);
}


void Viscoelastic::calcCt( const VectorXd & C)
{
	matHyper->calcCt(C);
	
	Ct = matHyper->Ct;
	Ctvol = matHyper->Ctvol;
	Ctiso = matHyper->Ctiso;
  
}

// ********************************************************************************************************

MatrixXd Viscoelastic::getTv(int i, const VectorXd & C, const MatrixXd & Gam)
{	return Tv.block(6*i,0,6,6);}

// ***

double Viscoelastic::getJ( const VectorXd & C)
{	return matHyper->getJ(C);}

//***

double Viscoelastic::getE()
{	return matHyper->getE();}

// ***

double Viscoelastic::getNu()
{	return matHyper->getNu();}

// ***

double Viscoelastic::getLambda()
{	return matHyper->getLambda();}

// ***

double Viscoelastic::getMu()
{	return matHyper->getMu();}

// ***

double Viscoelastic::getDensity()
{	return matHyper->density;}

// ***

VectorXd Viscoelastic::getTau()
{
	VectorXd Tau(nTau);
	
	for(int i=0;i<nTau;++i) Tau(i) = vTau[i];
	
	return Tau; 
}

// ***

VectorXd Viscoelastic::getBeta()
{
	VectorXd Beta(nTau);
	
	for(int i=0;i<nTau;++i) Beta(i) = vBeta[i];
	
	return Beta; 
}

// ***
VectorXd Viscoelastic::getVMu()
{       return vMu; }

// ***

VectorXd Viscoelastic::getEta()
{	return vEta;
}

// ***

int Viscoelastic::getNInternalVariable()
{	return nTau;}

// ***

int Viscoelastic::getMode()
{	return mode;}

// ***

int Viscoelastic::getRheological()
{       return rheological;}

// ***
